#define Analysis2_cxx
#include "Analysis2.h"
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TString.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <string>
#include <fstream>
#include <cmath>
#include "Escalas.h"
#include <iostream>

using namespace std;

void Analysis2::Loop()
{
  //gStyle->SetOptStat("nemruoisk");

  fChain->SetBranchStatus("*",1);
  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntriesFast();

  Long64_t nbytes = 0, nb = 0;

  const double GeV = 1000;
  const double TeV = 1000000;

  int nmetcf=13;
  int nbrpv=8;
  int s2dh=4; 

  // parameter cuts
  double srMUPT = 20*GeV;
  double srJJPT = 20*GeV;
  double srRJJM = 55*GeV;
  double srDRJJ = 2.9;
  double srDRJJMU = 1.1;
  double srTS = 0.25;
  double srMET = 130*GeV;
  double WMASS = 80.39*GeV;

  for(int i=0; i<nmetcf; i++){  
    metcf[i]=new TH1F(Form("metcf%d",i),Form("metcf%d",i),50,0,500);
  }


  for(int i=0; i<nbrpv; i++){
    if(i==4){
      met_brpv[i]=new TH1F(Form("met_brpv%d",i),Form("met_brpv%d",i),50,0,50);
      meff_brpv[i]=new TH1F(Form("meff_brpv%d",i),Form("meff_brpv%d",i),50,0,2);
      mt_brpv[i]=new TH1F(Form("mt_brpv%d",i),Form("mt_brpv%d",i),50,0,50);
      st_brpv[i]=new TH1F(Form("st_brpv%d",i),Form("st_brpv%d",i),50,0,1);
      ht_brpv[i]=new TH1F(Form("ht_brpv%d",i),Form("ht_brpv%d",i),50,0,2);
      met_over_meff_brpv[i]=new TH1F(Form("met_over_meff_brpv%d",i),Form("met_over_meff_brpv%d",i),50,0,1);
      mx_brpv[i]=new TH1F(Form("mx_brpv%d",i),Form("mx_brpv%d",i),100,0,500);
      mjj_brpv[i]=new TH1F(Form("mjj_brpv%d",i),Form("mjj_brpv%d",i),50,0,500);
      DRjj_brpv[i]=new TH1F(Form("DRjj_brpv%d",i),Form("DRjj_brpv%d",i),50,0,7);
      DRmujj_brpv[i]=new TH1F(Form("DRmujj_brpv%d",i),Form("DRmujj_brpv%d",i),50,0,7);
      mupt_brpv[i]=new TH1F(Form("mupt_brpv%d",i),Form("mupt_brpv%d",i),50,0,500);
      jjpt_brpv[i]=new TH1F(Form("jjpt_brpv%d",i),Form("jjpt_brpv%d",i),50,0,500);


      nm_brpv[i]=new TH1F(Form("nm_brpv%d",i),Form("nm_brpv%d",i),10,0,10);
      nj_brpv[i]=new TH1F(Form("nj_brpv%d",i),Form("nj_brpv%d",i),25,0,25);
      nw_brpv[i]=new TH1F(Form("nw_brpv%d",i),Form("nw_brpv%d",i),70,0,70);
      nx_brpv[i]=new TH1F(Form("nx_brpv%d",i),Form("nx_brpv%d",i),120,0,120);



      st_vs_ht_brpv[i]=new TH2F(Form("st_vs_ht_brpv%d",i),Form("st_vs_ht_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_met_brpv[i]=new TH2F(Form("st_vs_met_brpv%d",i),Form("st_vs_met_brpv%d",i),100/s2dh,0,50,100/s2dh,0,1);
      st_vs_meff_brpv[i]=new TH2F(Form("st_vs_meff_brpv%d",i),Form("st_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_mt_brpv[i]=new TH2F(Form("st_vs_mt_brpv%d",i),Form("st_vs_mt_brpv%d",i),100/s2dh,0,50,100/s2dh,0,1);

      mt_vs_met_brpv[i]=new TH2F(Form("mt_vs_met_brpv%d",i),Form("mt_vs_met_brpv%d",i),100/s2dh,0,50,100/s2dh,0,50);
      mt_vs_meff_brpv[i]=new TH2F(Form("mt_vs_meff_brpv%d",i),Form("mt_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,50);
      mt_vs_ht_brpv[i]=new TH2F(Form("mt_vs_ht_brpv%d",i),Form("mt_vs_ht_brpv%d",i),100/s2dh,0,2,100/s2dh,0,50);
      ht_vs_met_brpv[i]=new TH2F(Form("ht_vs_met_brpv%d",i),Form("ht_vs_met_brpv%d",i),100/s2dh,0,50,100/s2dh,0,2);


      st_vs_mt_vs_meff_brpv[i]=new TH3F(Form("st_vs_mt_vs_meff_brpv%d",i),Form("st_vs_mt_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,50,100/s2dh,0,1);
     
    }


    if(i==5 || i==6){
      met_brpv[i]=new TH1F(Form("met_brpv%d",i),Form("met_brpv%d",i),60,25,85);
      meff_brpv[i]=new TH1F(Form("meff_brpv%d",i),Form("meff_brpv%d",i),50,0,2);
      mt_brpv[i]=new TH1F(Form("mt_brpv%d",i),Form("mt_brpv%d",i),60,30,90);
      st_brpv[i]=new TH1F(Form("st_brpv%d",i),Form("st_brpv%d",i),50,0,1);
      ht_brpv[i]=new TH1F(Form("ht_brpv%d",i),Form("ht_brpv%d",i),50,0,2);
      met_over_meff_brpv[i]=new TH1F(Form("met_over_meff_brpv%d",i),Form("met_over_meff_brpv%d",i),50,0,1);
      mx_brpv[i]=new TH1F(Form("mx_brpv%d",i),Form("mx_brpv%d",i),100,0,500);
      mjj_brpv[i]=new TH1F(Form("mjj_brpv%d",i),Form("mjj_brpv%d",i),50,0,500);
      DRjj_brpv[i]=new TH1F(Form("DRjj_brpv%d",i),Form("DRjj_brpv%d",i),50,0,7);
      DRmujj_brpv[i]=new TH1F(Form("DRmujj_brpv%d",i),Form("DRmujj_brpv%d",i),50,0,7);
      mupt_brpv[i]=new TH1F(Form("mupt_brpv%d",i),Form("mupt_brpv%d",i),50,0,500);
      jjpt_brpv[i]=new TH1F(Form("jjpt_brpv%d",i),Form("jjpt_brpv%d",i),50,0,500);

      nm_brpv[i]=new TH1F(Form("nm_brpv%d",i),Form("nm_brpv%d",i),10,0,10);
      nj_brpv[i]=new TH1F(Form("nj_brpv%d",i),Form("nj_brpv%d",i),25,0,25);
      nw_brpv[i]=new TH1F(Form("nw_brpv%d",i),Form("nw_brpv%d",i),70,0,70);
      nx_brpv[i]=new TH1F(Form("nx_brpv%d",i),Form("nx_brpv%d",i),120,0,120);


      st_vs_ht_brpv[i]=new TH2F(Form("st_vs_ht_brpv%d",i),Form("st_vs_ht_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_met_brpv[i]=new TH2F(Form("st_vs_met_brpv%d",i),Form("st_vs_met_brpv%d",i),120/s2dh,24,84,100/s2dh,0,1);
      st_vs_meff_brpv[i]=new TH2F(Form("st_vs_meff_brpv%d",i),Form("st_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_mt_brpv[i]=new TH2F(Form("st_vs_mt_brpv%d",i),Form("st_vs_mt_brpv%d",i),120/s2dh,30,90,100/s2dh,0,1);

      mt_vs_met_brpv[i]=new TH2F(Form("mt_vs_met_brpv%d",i),Form("mt_vs_met_brpv%d",i),120/s2dh,24,84,120/s2dh,30,90);
      mt_vs_meff_brpv[i]=new TH2F(Form("mt_vs_meff_brpv%d",i),Form("mt_vs_meff_brpv%d",i),100/s2dh,0,2,120/s2dh,30,90);
      mt_vs_ht_brpv[i]=new TH2F(Form("mt_vs_ht_brpv%d",i),Form("mt_vs_ht_brpv%d",i),100/s2dh,0,2,120/s2dh,30,90);
      ht_vs_met_brpv[i]=new TH2F(Form("ht_vs_met_brpv%d",i),Form("ht_vs_met_brpv%d",i),120/s2dh,24,84,100/s2dh,0,2);


      st_vs_mt_vs_meff_brpv[i]=new TH3F(Form("st_vs_mt_vs_meff_brpv%d",i),Form("st_vs_mt_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,30,90,100/s2dh,0,1);
     
    }



    if(i!=4 && i!=5 && i!=6){
      met_brpv[i]=new TH1F(Form("met_brpv%d",i),Form("met_brpv%d",i),50,0,500);
      meff_brpv[i]=new TH1F(Form("meff_brpv%d",i),Form("meff_brpv%d",i),50,0,2);
      mt_brpv[i]=new TH1F(Form("mt_brpv%d",i),Form("mt_brpv%d",i),50,0,500);
      st_brpv[i]=new TH1F(Form("st_brpv%d",i),Form("st_brpv%d",i),50,0,1);
      ht_brpv[i]=new TH1F(Form("ht_brpv%d",i),Form("ht_brpv%d",i),50,0,2);
      met_over_meff_brpv[i]=new TH1F(Form("met_over_meff_brpv%d",i),Form("met_over_meff_brpv%d",i),50,0,1);
      mx_brpv[i]=new TH1F(Form("mx_brpv%d",i),Form("mx_brpv%d",i),100,0,500);
      mjj_brpv[i]=new TH1F(Form("mjj_brpv%d",i),Form("mjj_brpv%d",i),50,0,500);
      DRjj_brpv[i]=new TH1F(Form("DRjj_brpv%d",i),Form("DRjj_brpv%d",i),50,0,7);
      DRmujj_brpv[i]=new TH1F(Form("DRmujj_brpv%d",i),Form("DRmujj_brpv%d",i),50,0,7);
      mupt_brpv[i]=new TH1F(Form("mupt_brpv%d",i),Form("mupt_brpv%d",i),50,0,500);
      jjpt_brpv[i]=new TH1F(Form("jjpt_brpv%d",i),Form("jjpt_brpv%d",i),50,0,500);

      nm_brpv[i]=new TH1F(Form("nm_brpv%d",i),Form("nm_brpv%d",i),10,0,10);
      nj_brpv[i]=new TH1F(Form("nj_brpv%d",i),Form("nj_brpv%d",i),25,0,25);
      nw_brpv[i]=new TH1F(Form("nw_brpv%d",i),Form("nw_brpv%d",i),70,0,70);
      nx_brpv[i]=new TH1F(Form("nx_brpv%d",i),Form("nx_brpv%d",i),120,0,120);


      st_vs_ht_brpv[i]=new TH2F(Form("st_vs_ht_brpv%d",i),Form("st_vs_ht_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_met_brpv[i]=new TH2F(Form("st_vs_met_brpv%d",i),Form("st_vs_met_brpv%d",i),100/s2dh,0,500,100/s2dh,0,1);
      st_vs_meff_brpv[i]=new TH2F(Form("st_vs_meff_brpv%d",i),Form("st_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,1);
      st_vs_mt_brpv[i]=new TH2F(Form("st_vs_mt_brpv%d",i),Form("st_vs_mt_brpv%d",i),100/s2dh,0,500,100/s2dh,0,1);

      mt_vs_met_brpv[i]=new TH2F(Form("mt_vs_met_brpv%d",i),Form("mt_vs_met_brpv%d",i),100/s2dh,0,500,100/s2dh,0,500);
      mt_vs_meff_brpv[i]=new TH2F(Form("mt_vs_meff_brpv%d",i),Form("mt_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,500);
      mt_vs_ht_brpv[i]=new TH2F(Form("mt_vs_ht_brpv%d",i),Form("mt_vs_ht_brpv%d",i),100/s2dh,0,2,100/s2dh,0,500);
      ht_vs_met_brpv[i]=new TH2F(Form("ht_vs_met_brpv%d",i),Form("ht_vs_met_brpv%d",i),100/s2dh,0,500,100/s2dh,0,2);


      st_vs_mt_vs_meff_brpv[i]=new TH3F(Form("st_vs_mt_vs_meff_brpv%d",i),Form("st_vs_mt_vs_meff_brpv%d",i),100/s2dh,0,2,100/s2dh,0,500,100/s2dh,0,1);

    }
  }



  for(int i=0; i<nmetcf; i++){
    metcf[i]->Sumw2();
  }
  for(int i=0; i<nbrpv; i++){
    met_brpv[i]->Sumw2();
    meff_brpv[i]->Sumw2();
    mt_brpv[i]->Sumw2();
    st_brpv[i]->Sumw2();
    ht_brpv[i]->Sumw2();
    met_over_meff_brpv[i]->Sumw2();
    mx_brpv[i]->Sumw2();
    mjj_brpv[i]->Sumw2();
    DRjj_brpv[i]->Sumw2();
    DRmujj_brpv[i]->Sumw2();
    mupt_brpv[i]->Sumw2();
    jjpt_brpv[i]->Sumw2();
    nm_brpv[i]->Sumw2();
    nj_brpv[i]->Sumw2();
    nw_brpv[i]->Sumw2();
    nx_brpv[i]->Sumw2();


    st_vs_ht_brpv[i]->Sumw2();
    st_vs_met_brpv[i]->Sumw2();
    st_vs_meff_brpv[i]->Sumw2();
    st_vs_mt_brpv[i]->Sumw2();
    mt_vs_met_brpv[i]->Sumw2();
    mt_vs_meff_brpv[i]->Sumw2();
    mt_vs_ht_brpv[i]->Sumw2();
    ht_vs_met_brpv[i]->Sumw2();
    st_vs_mt_vs_meff_brpv[i]->Sumw2();
  }

  TFile* f = TFile::Open("plots.root","RECREATE");
  Escalas objesca; 
  double alfa;

  //**********************************************************************************************************************************
  //**********************************************************************************************************************************
  //**********************************************************************************************************************************
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    alfa = objesca.getEscala(RunNumber);   
    if(RunNumber==105009 || RunNumber==105010 || RunNumber==105011 || RunNumber==105012 ||
       RunNumber==105013 || RunNumber==105014 || RunNumber==105015 || RunNumber==105016) alfa*=0.55;

    //**************************
    metcf[0]->Fill(met/GeV,alfa);
    //**************************
    if(isGoodPV){
      metcf[1]->Fill(met/GeV,alfa);
      if(!wasBadJet){
        metcf[2]->Fill(met/GeV,alfa);
        if(MuN==1){
          metcf[3]->Fill(met/GeV,alfa);   
          if(MuPt->at(0)>20*GeV){
            metcf[4]->Fill(met/GeV,alfa);           
            if(JetN>=3){
              metcf[5]->Fill(met/GeV,alfa);
              if(JetPt->at(0)>60*GeV){
                metcf[6]->Fill(met/GeV,alfa);
                if(JetPt->at(1)>30*GeV){
                  metcf[7]->Fill(met/GeV,alfa);
                  if(JetPt->at(2)>30*GeV){
                    metcf[8]->Fill(met/GeV,alfa);
                    if(DeltaPhiMetJet->at(0)>0.2 && DeltaPhiMetJet->at(1)>0.2 && DeltaPhiMetJet->at(2)>0.2){


                      metcf[9]->Fill(met/GeV,alfa); 
                      if(Muz0_exPV->at(0)<10.0){
                        // At this point I believe that this is the KS of the 1l3j analysis.
                        // The Event veto for crack electrons requested by the twiki page, is also changed by a removal electron.
                        metcf[10]->Fill(met/GeV,alfa);
                      
                        if(!wasjptlet30){
                          metcf[11]->Fill(met/GeV,alfa);
                          if(MuTrigger){
                            metcf[12]->Fill(met/GeV,alfa);
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }

    // Here start the thesis analysis
    if(isGoodPV){
      if(!wasBadJet){
        if(MuN>=1){
          if(JetN>=3){
            if(JetPt->at(0)>60*GeV){
              if(JetPt->at(1)>30*GeV){
                if(JetPt->at(2)>30*GeV){
                  if(DeltaPhiMetJet->at(0)>0.2 && DeltaPhiMetJet->at(1)>0.2 && DeltaPhiMetJet->at(2)>0.2){
                    if(!wasCosmicMuon){
                      nbjets=0;
                      for(unsigned int i=0; i<3; i++){
                        if(JetFlavorWeightSV0->at(i)>5.72) nbjets++; 
                      }

                      met_brpv[0]->Fill(met/GeV,alfa);
                      meff_brpv[0]->Fill(em_muonjetjetjet/TeV,alfa);
                      st_brpv[0]->Fill(ts,alfa);
                      ht_brpv[0]->Fill(ht_muonjetjetjet/TeV,alfa);
                      met_over_meff_brpv[0]->Fill(met/em_muonjetjetjet,alfa);
                      st_vs_ht_brpv[0]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
                      st_vs_met_brpv[0]->Fill(met/GeV,ts,alfa);
                      st_vs_meff_brpv[0]->Fill(em_muonjetjetjet/TeV,ts,alfa);
                      ht_vs_met_brpv[0]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[0]->Fill(MuN,alfa);
                      nj_brpv[0]->Fill(JetN,alfa);
                      nw_brpv[0]->Fill(m_jj->size(),alfa);
                      nx_brpv[0]->Fill(m_Wmu->size(),alfa);

                      
                      for(int i=0; i<MuN; i++){
			mupt_brpv[0]->Fill(MuPt->at(i)/GeV,alfa);
                        mt_brpv[0]->Fill(MTs->at(i)/GeV,alfa);
                        mt_vs_met_brpv[0]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
                        mt_vs_meff_brpv[0]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                        mt_vs_ht_brpv[0]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                        st_vs_mt_brpv[0]->Fill(MTs->at(i)/GeV,ts,alfa);
                        st_vs_mt_vs_meff_brpv[0]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                      }
                      for(unsigned int i=0; i<m_Wmu->size(); i++){
                        mx_brpv[0]->Fill(m_Wmu->at(i)/GeV,alfa);
                        DRmujj_brpv[0]->Fill(DeltaR_Wmu->at(i),alfa);
                      }
                      for(unsigned int i=0; i<m_jj->size(); i++){
                        mjj_brpv[0]->Fill(m_jj->at(i)/GeV,alfa);
                        DRjj_brpv[0]->Fill(DeltaR_jj->at(i),alfa);
                        jjpt_brpv[0]->Fill(pt_jj->at(i)/GeV,alfa); 
                      } 




                      if(em_muonjetjetjet>500*GeV){
                        met_brpv[1]->Fill(met/GeV,alfa);
                        meff_brpv[1]->Fill(em_muonjetjetjet/TeV,alfa);
                        st_brpv[1]->Fill(ts,alfa);
                        ht_brpv[1]->Fill(ht_muonjetjetjet/TeV,alfa);
                        met_over_meff_brpv[1]->Fill(met/em_muonjetjetjet,alfa);
                        st_vs_ht_brpv[1]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
                        st_vs_met_brpv[1]->Fill(met/GeV,ts,alfa);
                        st_vs_meff_brpv[1]->Fill(em_muonjetjetjet/TeV,ts,alfa);
                        ht_vs_met_brpv[1]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                        nm_brpv[1]->Fill(MuN,alfa);
                        nj_brpv[1]->Fill(JetN,alfa);
                        nw_brpv[1]->Fill(m_jj->size(),alfa);
                        nx_brpv[1]->Fill(m_Wmu->size(),alfa);

                        for(int i=0; i<MuN; i++){
      			  mupt_brpv[1]->Fill(MuPt->at(i)/GeV,alfa);
                          mt_brpv[1]->Fill(MTs->at(i)/GeV,alfa);
                          mt_vs_met_brpv[1]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
                          mt_vs_meff_brpv[1]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                          mt_vs_ht_brpv[1]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                          st_vs_mt_brpv[1]->Fill(MTs->at(i)/GeV,ts,alfa);
                          st_vs_mt_vs_meff_brpv[1]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                        }
                        for(unsigned int i=0; i<m_Wmu->size(); i++){
                          mx_brpv[1]->Fill(m_Wmu->at(i)/GeV,alfa);
                          DRmujj_brpv[1]->Fill(DeltaR_Wmu->at(i),alfa);
                        }
                        for(unsigned int i=0; i<m_jj->size(); i++){
                          mjj_brpv[1]->Fill(m_jj->at(i)/GeV,alfa);
                          DRjj_brpv[1]->Fill(DeltaR_jj->at(i),alfa);
                          jjpt_brpv[1]->Fill(pt_jj->at(i)/GeV,alfa); 
                        } 

			if(met/em_muonjetjetjet>0.25){
			  met_brpv[2]->Fill(met/GeV,alfa);
			  meff_brpv[2]->Fill(em_muonjetjetjet/TeV,alfa);
			  st_brpv[2]->Fill(ts,alfa);
			  ht_brpv[2]->Fill(ht_muonjetjetjet/TeV,alfa);
			  met_over_meff_brpv[2]->Fill(met/em_muonjetjetjet,alfa);
			  st_vs_ht_brpv[2]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
			  st_vs_met_brpv[2]->Fill(met/GeV,ts,alfa);
			  st_vs_meff_brpv[2]->Fill(em_muonjetjetjet/TeV,ts,alfa);
			  ht_vs_met_brpv[2]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[2]->Fill(MuN,alfa);
                      nj_brpv[2]->Fill(JetN,alfa);
                      nw_brpv[2]->Fill(m_jj->size(),alfa);
                      nx_brpv[2]->Fill(m_Wmu->size(),alfa);

			  for(int i=0; i<MuN; i++){
                            mupt_brpv[2]->Fill(MuPt->at(i)/GeV,alfa);
			    mt_brpv[2]->Fill(MTs->at(i)/GeV,alfa);
			    mt_vs_met_brpv[2]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
			    mt_vs_meff_brpv[2]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
			    mt_vs_ht_brpv[2]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
			    st_vs_mt_brpv[2]->Fill(MTs->at(i)/GeV,ts,alfa);
			    st_vs_mt_vs_meff_brpv[2]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
			  }
                          for(unsigned int i=0; i<m_Wmu->size(); i++){
                            mx_brpv[2]->Fill(m_Wmu->at(i)/GeV,alfa);
                            DRmujj_brpv[2]->Fill(DeltaR_Wmu->at(i),alfa);
                          }
                          for(unsigned int i=0; i<m_jj->size(); i++){
                            mjj_brpv[2]->Fill(m_jj->at(i)/GeV,alfa);
                            DRjj_brpv[2]->Fill(DeltaR_jj->at(i),alfa);
                            jjpt_brpv[2]->Fill(pt_jj->at(i)/GeV,alfa); 
                          } 


			} // met/em_muonjetjejjet>0.25
		      }// em_muonjetjetjet>500*GeV

                      // signal region 
  
                      if( ((ts>0.125) && (ts>1.0-2.1875*(em_muonjetjetjet/TeV-0.4))) || ((ts<=0.125) && (em_muonjetjetjet/TeV>0.8)) ){
                        met_brpv[3]->Fill(met/GeV,alfa);
                        meff_brpv[3]->Fill(em_muonjetjetjet/TeV,alfa);
                        st_brpv[3]->Fill(ts,alfa);
                        ht_brpv[3]->Fill(ht_muonjetjetjet/TeV,alfa);
                        met_over_meff_brpv[3]->Fill(met/em_muonjetjetjet,alfa);
                        st_vs_ht_brpv[3]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
                        st_vs_met_brpv[3]->Fill(met/GeV,ts,alfa);
                        st_vs_meff_brpv[3]->Fill(em_muonjetjetjet/TeV,ts,alfa);
                        ht_vs_met_brpv[3]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[3]->Fill(MuN,alfa);
                      nj_brpv[3]->Fill(JetN,alfa);
                      nw_brpv[3]->Fill(m_jj->size(),alfa);
                      nx_brpv[3]->Fill(m_Wmu->size(),alfa);

                        for(int i=0; i<MuN; i++){
                          mupt_brpv[3]->Fill(MuPt->at(i)/GeV,alfa);
                          mt_brpv[3]->Fill(MTs->at(i)/GeV,alfa);
                          mt_vs_met_brpv[3]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
                          mt_vs_meff_brpv[3]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                          mt_vs_ht_brpv[3]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                          st_vs_mt_brpv[3]->Fill(MTs->at(i)/GeV,ts,alfa);
                          st_vs_mt_vs_meff_brpv[3]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                        }
                        for(unsigned int i=0; i<m_jj->size(); i++){
			  mjj_brpv[3]->Fill(m_jj->at(i)/GeV,alfa);
		          DRjj_brpv[3]->Fill(DeltaR_jj->at(i),alfa);
		          jjpt_brpv[3]->Fill(pt_jj->at(i)/GeV,alfa);
                        } 
                        for(unsigned int i=0; i<m_Wmu->size(); i++){
		          mx_brpv[3]->Fill(m_Wmu->at(i)/GeV,alfa);
			  DRmujj_brpv[3]->Fill(DeltaR_Wmu->at(i),alfa);
                        }











                        if(DEBUG) cout << "before v_1 fills" <<endl;
                        for(unsigned int i=0; i<m_Wmu->size(); i++){
                          unsigned int j=i/m_jj->size(); // muon index
                          unsigned int k=i/MuN; // jj index
			    if(MuPt->at(j)>srMUPT){
			      if(fabs(m_jj->at(k)-WMASS)<srRJJM){
				if(pt_jj->at(k)>srJJPT){
				  if(DeltaR_jj->at(k)<srDRJJ){
				    if(DeltaR_Wmu->at(i)<srDRJJMU){
                                      if(ts>srTS){
                                        if(met>srMET){
				          met_brpv[7]->Fill(met/GeV,alfa);
				          meff_brpv[7]->Fill(em_muonjetjetjet/TeV,alfa);
				          st_brpv[7]->Fill(ts,alfa);
				          ht_brpv[7]->Fill(ht_muonjetjetjet/TeV,alfa);
				          met_over_meff_brpv[7]->Fill(met/em_muonjetjetjet,alfa);

				          st_vs_ht_brpv[7]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
				          st_vs_met_brpv[7]->Fill(met/GeV,ts,alfa);
				          st_vs_meff_brpv[7]->Fill(em_muonjetjetjet/TeV,ts,alfa);
				          ht_vs_met_brpv[7]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[7]->Fill(MuN,alfa);
                      nj_brpv[7]->Fill(JetN,alfa);
                      nw_brpv[7]->Fill(m_jj->size(),alfa);
                      nx_brpv[7]->Fill(m_Wmu->size(),alfa);

				          break;
                                        }
                                      }
				    }
				  }
				}
			      }
			    }
                        } // endfor
                        if(DEBUG) cout << "before v_mu fills" <<endl;
                        for(int j=0; j<MuN; j++){
                          if(MuPt->at(j)>srMUPT){
                            for(unsigned int i=0; i<m_jj->size(); i++){
                              if(fabs(m_jj->at(i)-WMASS)<srRJJM){
                                if(pt_jj->at(i)>srJJPT){
                                  if(DeltaR_jj->at(i)<srDRJJ){
                                    if( DeltaR_Wmu->at(i+j*m_jj->size())<srDRJJMU ){
                                      if(ts>srTS){
                                        if(met>srMET){
					  mupt_brpv[7]->Fill(MuPt->at(j)/GeV,alfa);
					  mt_brpv[7]->Fill(MTs->at(j)/GeV,alfa);

					  mt_vs_met_brpv[7]->Fill(met/GeV,MTs->at(j)/GeV,alfa);
					  mt_vs_meff_brpv[7]->Fill(em_muonjetjetjet/TeV,MTs->at(j)/GeV,alfa);
					  mt_vs_ht_brpv[7]->Fill(ht_muonjetjetjet/TeV,MTs->at(j)/GeV,alfa);
					  st_vs_mt_brpv[7]->Fill(MTs->at(j)/GeV,ts,alfa);
					  st_vs_mt_vs_meff_brpv[7]->Fill(em_muonjetjetjet/TeV,MTs->at(j)/GeV,ts,alfa);
                                        }
                                      }
                                    }
                                  }
                                }
                              }
                            }
                          }
                        }
                        for(unsigned int i=0; i<m_jj->size(); i++){
                          if(fabs(m_jj->at(i)-WMASS)<srRJJM){
                            if(pt_jj->at(i)>srJJPT){
			      if(DeltaR_jj->at(i)<srDRJJ){ 
			        for(int j=0; j<MuN; j++){
                                  if(MuPt->at(j)>srMUPT){
				    if(DeltaR_Wmu->at(i+j*m_jj->size())<srDRJJMU){
                                      if(ts>srTS){
                                        if(met>srMET){
					  mjj_brpv[7]->Fill(m_jj->at(i)/GeV,alfa);
					  DRjj_brpv[7]->Fill(DeltaR_jj->at(i),alfa);
					  jjpt_brpv[7]->Fill(pt_jj->at(i)/GeV,alfa);
					  break; // to avoid degeneracy on n_{muon}
                                        }
                                      }
			 	    }
                                  }
			        }
                              } 
			    }
                          }
                        } 
                        for(unsigned int i=0; i<m_Wmu->size(); i++){
                          unsigned int j=i/m_jj->size(); // muon index
                          unsigned int k=i/MuN; // jj index
                          if(MuPt->at(j)>srMUPT){
			    if(fabs(m_jj->at(k)-WMASS)<srRJJM){
                              if(pt_jj->at(k)>srJJPT){
			        if(DeltaR_jj->at(k)<srDRJJ){
				  if(DeltaR_Wmu->at(i)<srDRJJMU){
                                    if(ts>srTS){
                                      if(met>srMET){
				        mx_brpv[7]->Fill(m_Wmu->at(i)/GeV,alfa);
				        DRmujj_brpv[7]->Fill(DeltaR_Wmu->at(i),alfa);
                                      }
                                    }
				  }
			        }
                              }
			    }
                          }
                        }











                      } // cut in the meff,ts plane

                      // QCD CR
                      unsigned int n_pev_mtlt40=0;
                      for(unsigned int j=0; j<MTs->size(); j++){
                        if(MTs->at(j)<40*GeV) n_pev_mtlt40++;
                      }
                      if(met<40*GeV && n_pev_mtlt40>0){                    
                        met_brpv[4]->Fill(met/GeV,alfa);
                        meff_brpv[4]->Fill(em_muonjetjetjet/TeV,alfa);
                        st_brpv[4]->Fill(ts,alfa);
                        ht_brpv[4]->Fill(ht_muonjetjetjet/TeV,alfa);
                        met_over_meff_brpv[4]->Fill(met/em_muonjetjetjet,alfa);
                        st_vs_ht_brpv[4]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
                        st_vs_met_brpv[4]->Fill(met/GeV,ts,alfa);
                        st_vs_meff_brpv[4]->Fill(em_muonjetjetjet/TeV,ts,alfa);
                        ht_vs_met_brpv[4]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[4]->Fill(MuN,alfa);
                      nj_brpv[4]->Fill(JetN,alfa);
                      nw_brpv[4]->Fill(m_jj->size(),alfa);
                      nx_brpv[4]->Fill(m_Wmu->size(),alfa);

                        for(unsigned int i=0; i<m_jj->size(); i++){
                          mjj_brpv[4]->Fill(m_jj->at(i)/GeV,alfa);
                          DRjj_brpv[4]->Fill(DeltaR_jj->at(i),alfa);
                          jjpt_brpv[4]->Fill(pt_jj->at(i)/GeV,alfa); 
                        } 
                        for(int i=0; i<MuN; i++){
                          if(MTs->at(i)<40*GeV){
                            mupt_brpv[4]->Fill(MuPt->at(i)/GeV,alfa);
                            mt_brpv[4]->Fill(MTs->at(i)/GeV,alfa);
                            mt_vs_met_brpv[4]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
                            mt_vs_meff_brpv[4]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                            mt_vs_ht_brpv[4]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                            st_vs_mt_brpv[4]->Fill(MTs->at(i)/GeV,ts,alfa);
                            st_vs_mt_vs_meff_brpv[4]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                          }
                        }
                        for(unsigned int i=0; i<m_Wmu->size(); i++){
                          unsigned int j=i/m_jj->size();
                          if(MTs->at(j)<40*GeV){
                            mx_brpv[4]->Fill(m_Wmu->at(i)/GeV,alfa);
                            DRmujj_brpv[4]->Fill(DeltaR_Wmu->at(i),alfa);
                          }
                        }


                      }
                      // W/T CR
                      int n_pev_mtge40lt80=0;
                      for(unsigned int j=0; j<MTs->size(); j++){
                        if(MTs->at(j)>=40*GeV && MTs->at(j)<80*GeV) n_pev_mtge40lt80++;
                      }

                      if(met>=30*GeV && met<80*GeV && n_pev_mtge40lt80>0){
                        if(nbjets==0){
			  met_brpv[5]->Fill(met/GeV,alfa);
			  meff_brpv[5]->Fill(em_muonjetjetjet/TeV,alfa);
			  st_brpv[5]->Fill(ts,alfa);
			  ht_brpv[5]->Fill(ht_muonjetjetjet/TeV,alfa);
			  met_over_meff_brpv[5]->Fill(met/em_muonjetjetjet,alfa);
			  st_vs_ht_brpv[5]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
			  st_vs_met_brpv[5]->Fill(met/GeV,ts,alfa);
			  st_vs_meff_brpv[5]->Fill(em_muonjetjetjet/TeV,ts,alfa);
			  ht_vs_met_brpv[5]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[5]->Fill(MuN,alfa);
                      nj_brpv[5]->Fill(JetN,alfa);
                      nw_brpv[5]->Fill(m_jj->size(),alfa);
                      nx_brpv[5]->Fill(m_Wmu->size(),alfa);

			  for(unsigned int i=0; i<m_jj->size(); i++){
			    mjj_brpv[5]->Fill(m_jj->at(i)/GeV,alfa);
			    DRjj_brpv[5]->Fill(DeltaR_jj->at(i),alfa);
			    jjpt_brpv[5]->Fill(pt_jj->at(i)/GeV,alfa); 
			  } 
			  for(int i=0; i<MuN; i++){
                            if(MTs->at(i)>=40*GeV && MTs->at(i)<80*GeV){
                              mupt_brpv[5]->Fill(MuPt->at(i)/GeV,alfa);
			      mt_brpv[5]->Fill(MTs->at(i)/GeV,alfa);
			      mt_vs_met_brpv[5]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
			      mt_vs_meff_brpv[5]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
			      mt_vs_ht_brpv[5]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
			      st_vs_mt_brpv[5]->Fill(MTs->at(i)/GeV,ts,alfa);
			      st_vs_mt_vs_meff_brpv[5]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                            }
			  }
			  for(unsigned int i=0; i<m_Wmu->size(); i++){
                            unsigned int j=i/m_jj->size();
                            if(MTs->at(j)>=40*GeV && MTs->at(j)<80*GeV){
			      mx_brpv[5]->Fill(m_Wmu->at(i)/GeV,alfa);
			      DRmujj_brpv[5]->Fill(DeltaR_Wmu->at(i),alfa);
                            }
			  }
                        }
                        if(nbjets>0){
                          met_brpv[6]->Fill(met/GeV,alfa);
                          meff_brpv[6]->Fill(em_muonjetjetjet/TeV,alfa);
                          st_brpv[6]->Fill(ts,alfa);
                          ht_brpv[6]->Fill(ht_muonjetjetjet/TeV,alfa);
                          met_over_meff_brpv[6]->Fill(met/em_muonjetjetjet,alfa);
                          st_vs_ht_brpv[6]->Fill(ht_muonjetjetjet/TeV,ts,alfa);
                          st_vs_met_brpv[6]->Fill(met/GeV,ts,alfa);
                          st_vs_meff_brpv[6]->Fill(em_muonjetjetjet/TeV,ts,alfa);
                          ht_vs_met_brpv[6]->Fill(met/GeV,ht_muonjetjetjet/TeV,alfa);
                      nm_brpv[6]->Fill(MuN,alfa);
                      nj_brpv[6]->Fill(JetN,alfa);
                      nw_brpv[6]->Fill(m_jj->size(),alfa);
                      nx_brpv[6]->Fill(m_Wmu->size(),alfa);
                          for(unsigned int i=0; i<m_jj->size(); i++){
                            mjj_brpv[6]->Fill(m_jj->at(i)/GeV,alfa);
                            DRjj_brpv[6]->Fill(DeltaR_jj->at(i),alfa);
                            jjpt_brpv[6]->Fill(pt_jj->at(i)/GeV,alfa); 
                          } 
                          for(int i=0; i<MuN; i++){
                            if(MTs->at(i)>=40*GeV && MTs->at(i)<80*GeV){
                              mupt_brpv[6]->Fill(MuPt->at(i)/GeV,alfa);
                              mt_brpv[6]->Fill(MTs->at(i)/GeV,alfa);
                              mt_vs_met_brpv[6]->Fill(met/GeV,MTs->at(i)/GeV,alfa);
                              mt_vs_meff_brpv[6]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                              mt_vs_ht_brpv[6]->Fill(ht_muonjetjetjet/TeV,MTs->at(i)/GeV,alfa);
                              st_vs_mt_brpv[6]->Fill(MTs->at(i)/GeV,ts,alfa);
                              st_vs_mt_vs_meff_brpv[6]->Fill(em_muonjetjetjet/TeV,MTs->at(i)/GeV,ts,alfa);
                            }
                          }
                          for(unsigned int i=0; i<m_Wmu->size(); i++){
                            unsigned int j=i/m_jj->size();
                            if(MTs->at(j)>=40*GeV && MTs->at(j)<80*GeV){
                              mx_brpv[6]->Fill(m_Wmu->at(i)/GeV,alfa);
                              DRmujj_brpv[6]->Fill(DeltaR_Wmu->at(i),alfa);
                            }
                          }


                        }

                      } // W/T CR




		    } // !wasCosmicMuon
		  }
		}
	      }
	    }
	  }
	}
      }
    }


  }
  //**********************************************************************************************************************************
  //**********************************************************************************************************************************
  //**********************************************************************************************************************************

  f->cd();
  for(int i=0; i<nmetcf; i++){
    metcf[i]->Write();
  }
  for(int i=0; i<nbrpv; i++){

    met_brpv[i]->Write();
    meff_brpv[i]->Write();
    mt_brpv[i]->Write();
    st_brpv[i]->Write();
    ht_brpv[i]->Write();
    met_over_meff_brpv[i]->Write();
    mx_brpv[i]->Write();
    mjj_brpv[i]->Write();
    DRjj_brpv[i]->Write();
    DRmujj_brpv[i]->Write();
    mupt_brpv[i]->Write();
    jjpt_brpv[i]->Write();
    nm_brpv[i]->Write();
    nj_brpv[i]->Write();
    nw_brpv[i]->Write();
    nx_brpv[i]->Write();



    st_vs_ht_brpv[i]->Write();
    st_vs_met_brpv[i]->Write();
    st_vs_meff_brpv[i]->Write();
    st_vs_mt_brpv[i]->Write();
    mt_vs_met_brpv[i]->Write();
    mt_vs_meff_brpv[i]->Write();
    mt_vs_ht_brpv[i]->Write();
    ht_vs_met_brpv[i]->Write();
    st_vs_mt_vs_meff_brpv[i]->Write();
  }

  f->Close();

}

Analysis2::Analysis2( string inputFile )
{

  TChain* ch = new TChain("Nt");
  string fileName;
  ifstream ffin( inputFile.data() );
  while( ffin ){
    getline( ffin , fileName );
    if( fileName.length() > 0 ) ch -> AddFile( fileName.data() );
  }
  Init(ch);
  DEBUG = false;



}


